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CS| ' We study the convective instability of the melt during the initial transient in a di- 

Q \ rectional solidification experiment in a vertical configuration. We obtain analytically 

^ ■ the dispersion relation, and perform an additional asymptotic expansion for large 

solutal Rayleigh number that permits a simpler analytical analysis and a better nu- 
merical behavior. We find a transient instability, i.e. a regime in which the system 
j-^' destabilizes during the transient whereas the final unperturbed steady state is stable. 

^ ! This could be relevant to growth mode predictions in solidification. 
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I. INTRODUCTION 

When a melt solidifies the advancing solid-liquid interface can develop very complex 
patterns, which induce inhomogeneities in the solid phase. The resulting microstructure is 
largely responsible for the final properties of the solid. ^ A significative number of methods 
employed in Materials Science involve solidification of a melt, and hence prediction of the 
detailed dynamics of the solidification front turns out to be of capital importance. Most 
solidification processes of applied interest can be modelized by the directional solidification 
configuration. In it, a mixture or alloy is pulled at a controlled velocity inside a thermal 
gradient. This setup permits to control to a certain degree the growth velocity of the solid 
phase, which becomes a control parameter of the problem.- Even ignoring convection in 
the melt (or supressing it in experiments) the phenomenology obtained with this setup 
is very rich, with growth of cells and dendrites (with the presence of sidebranching) , but 
also with less usual structures such as seaweeds or doublons, and non steady behaviors 
such as traveling waves, banding, etc. One interesting point is that these processes can be 
history dependent. For instance, in dendritic growth, solutions with different wavelengths 
are possible for a given growth condition,- and the final one is dynamically selected during 
the initial transient states of the solidification process.'^'^ The fact is that after setting the 
final value of the pulling velocity, there is a transient in which a layer of solute is built 
up ahead of the advancing front, during which the instability can occur .-i^ In the current 
understanding of the problem this transient determines the result of the selection process.-i^ 
Indeed, calculations performed by using the instability of the steady state are known to yield 
wrong results on the selected wavelength when comparing to experiments.-'^ 

While thermal and concentration gradients can induce convection in the melt, this effect 
has received little attention in time- dependent solidification problems, but it has nevertheless 
been studied by using different approximations in some solidification setup3i2,"i^. 

Convection is unavoidable in the case of lateral heating, i.e. when the thermal gradient is 
placed horizontally. On the contrary we will center here on a vertical configuration, in which 
the solidifying sample is immersed in a vertical thermal gradient and is slowly extracted 
from below. In this case thermal gradient is stabilizing, but accumulation of solute ahead 
the solidification front, depending on its relative density, can induce an instability of the 
Rayleigh-Benard type. For this later configuration there have been numerous works studying 



the convective instability of the steady state. Critical parameters and wavelength of the 
instability were found by Hurle et alM A weakly nonlinear analysis by Jenkins^ derived 
an amplitude equation that supported hexagonal solutions. Riley and Davis^^ derived an 
asymptotic scaling law for the critical wavenumber in the same system. Also, a weakly 
non-linear analysis was performed resulting in a Sivashinsky equation. Relatively complete 
studies of the problem of convection in the melt was undertaken by Impey et al- for 
steady solutions and LeMarec et alM for time-dependent states, finding chaotic dynamics. 
LeMarec also undertook similar simulations in three dimensions,— focusing in the prediction 
of actual macrosegregation patterns in the final crystal. Effects of the solidification transient 
on the convective instability were identified by the experiments of Jamgotchian et alM 
who obtained segregation patterns induced by thermosolutal convection during a directional 
solidification experiment. Here the size of the convective cells did not agree with the steady- 
state wavelength for the hydrodynamic instability, but was found to be consistent with 
estimations obtained taking into account the transient dynamics. 



In this paper, we address the question of the first convective instability of the melt during 
the initial transient of a solidification experiment in the vertical configuration. To this end, 
it is necessary to employ a time dependent solution of the solidification problem, i.e. the 
free boundary dynamics of a planar front with the appropriate boundary conditions at the 
interface and the diffusion equation for the solute in the melt. We will use the approximate 
solution by Warren and Langer,- which is known to be remarkably good when compared to 
the exact solution.— The analytical structure of such solution adapts well for the stability 
analysis of the flow, which will be performed semi-analytically as a function of time during 
all the transient. In doing so, time is treated as a parameter, thus performing a quasi-static 
or adiabatic approximation^ii^. Furthermore, an asymptotic expansion for large Rayleigh 
number will permit to gain physical insight on the obtained solutions. One interesting 
possibility will be the appearance of a transient instability, i.e. an instability occurring 
during the transient stage but which is not present in the steady state. Due to the fact 
that any present flow will induce a deformation of the solidification front, such transient 
instability is also expected to affect the final selected growth mode. We will find a regime 
in which such instability appears. 




(a) 



(b) 




FIG. 1. (a) Setup of the system. The sample is slowly pulled from a furnace, represented as a 
thermal gradient, (b) Equilibrium configuration and its corresponding solute profile, for i^ < 1 (c) 
Steady state configuration and concentration profile. The possibility of convection in the melt is 
also sketched. 



II. THE MODEL 



In the vertical configuration of a directional solidification experiment we place a sample 
with a mixture or alloy into a vertical constant thermal gradient as shown in Fig. [TK (it 
is assumed to be time independent owing to the large thermal conductivity. This is the 
frozen temperature approximation^^^) . The sample is initially at rest with a constant solute 
concentration C^o in the molten phase (Fig. Wp)- The interface is planar and its initial 
position is given by the equilibrium temperature at that concentration following the phase 
diagram of the mixture-. The thermal gradient is stabilizing and consequently the melt is 
quiescent. At time t = we start to pull from the solid at a constant speed Vp and, if 
the interface is stable, a solutal boundary layer is built ahead of the interface in the steady 
state (Fig. [lb). Note the recoil of the solidification front, due to the change in the interface 
concentration during the process, and the fact that this means a change in the coexistence 
temperature. 



A. Basic equations 

Mathematical description of this system will consist in the Navier-Stokes equations for the 
fluid flow, with the appropriate boundary conditions, in the the Boussinesq approximation,— 
i.e. we will assume the density as a constant mean value except in the buoyancy term. 
Furthermore we will consider buoyancy effects coming from density changes due to solute 
concentration only. Effect of thermal buoyancy is stabilizing in the present configuration 
and would only be of importance at small velocities and solute concentrations.^^ For the 
solidification process we will use the so-called minimal model of directional solidification,- 
which consists in the diffusion equation for the concentration, supplemented with the bound- 
ary conditions (solute conservation and local equilibrium) at the (free moving) interface. All 
equations will be written exclusively in the liquid phase, assuming that neither flow or 
diffusion are possible in the solid phase. 

The change in density will depend on concentration through a linear approximation: 

p = P0{l + a{C-Cref)), (1) 

where po is a reference density taken at the value of the reference concentration Cref, and a 
is the solute expansion coefficient, which will be assumed to be negative, corresponding to 
a buoyant solute. 

With the previous approximations taken into account, the governing equations for the 
velocity field v and the solute concentration in the bulk C will take the following form, 
written in the gradient frame: 

dtv - t;p 9,v + (v ■ V) V = g (1 + a (C - Cref)) - Vvr + z/VV, (2) 

Vv = 0, (3) 

dtC-Vpd,C = DV^C, (4) 

where Vp is the value of the pulling speed, g is the acceleration of gravity, D is the solute 
diffusion coefficient in the liquid phase, u is the kinematic viscosity and vr is the pressure 
reduced by the reference density. Axis z is placed in the direction of the thermal gradient. 
The previous equations are nondimensionalized by using the following scalings (tildes 



denote non-dimensional parameters): 

r = Idt, t = TDt, V = fpV, 



(5) 



C = CooC, IT = VpTT. 



Here Coo is the initial concentration of the melt. We also use Id = D/vp and td = D/Vp 
as the diffusion length and time, respectively. Dropping tildes and assuming that we have 
eliminated the constant terms by a suitable redefinition of the pressure, taking twice the 
curl of Eq. [2] and after some manipulation, the following equations are obtained: 

9tV2v-9,V2v + V2((v- V)v)-VN(v) = ScV^V\ 

+ R^s-^^ [idle + die) e, - d^dyCey - d^d^Ce^) , (6) 

dtC - d,c = v'c, (7) 

where 

N(v) = {d.v,){d,v,) (8) 

(summation over both indices). The z component of Eq. |6]can be rewritten as: 

dtV\, - d-y%, + V^ ((v ■ V) V,) - 9,N(v) = ScV^V%-, (9) 

+ Ras^^{dlC + dlc). 

Rayleigh number Ras and Schmidt number Sc are defined as:— 

Sc = ^. (11) 



Here we have introduced the segregation constant K. Assuming a dilute alloy the phase 
diagram is taken as linear, and the concentrations of the coexistent phases at the interface 
are related through the segregation constant, i.e. CsoUd = KC at the interface, which we 
assume to be i^ < 1. Note also that Ras is positive since we have a buoyant solute (a < 0). 
Pressure can be recovered from the solution of Eq. [H] by using the following equation: 

V^TT = Ras-^d,C + N(v). (12) 

i — A 



B. Boundary conditions 

Boundary conditions have to be established for both velocity and concentration fields. 
The velocity field has no-slip boundary conditions at the solidifying interface. Furthermore 
its normal velocity is zero in the frame moving with the sample. If velocity (in the gra- 
dient frame) is redefined subtracting the pulling velocity, it can be seen that all velocity 
components and their first spatial derivatives do vanish at the interface: 

v(x,2/,C(t)) = 0, (13) 

^^^.■L=CW = 0' Vz,j = 1...3, (14) 

being ({t) the z coordinate of the planar interface at time t. Furthermore, since we neglect 
density changes during the phase transition, the redefined velocity will also vanish at infinity, 
that is 

V = at z = CO. (15) 

If pressure has to be recovered from Eq. [121 a condition for the normal derivative of the 
pressure at the interface can be easily found by taking the z component of Eq. |2l 

For the concentration field there are two boundary conditions at the interface. First of 
all, we have the solute conservation equation, which in dimensional variables is written as 

Dd,C=iK-l)Cat). (16) 

The first term stands for the diffusion fiux and the term in the right accounts for the solute 
expelled by the front due to the different equilibrium concentrations of both phases, ({t) is 
the interface velocity. We have implicitly assumed that there is no solute diffusivity in the 
solid. The other boundary condition is the Gibbs- Thompson equation, which corresponds to 
the hypothesis of local equilibrium at the interface and for a planar front reads^ 

Ti = TM + mLC, (17) 

where Tj is the interface temperature, Tm is the melting temperature of the pure solvent, and 
rriL is the (negative) coexistence liquidus slope in the linearized equilibrium phase diagram of 
the alloy. This last equation will permit to relate the concentration to the interface position 
in the gradient frame. We take the origin z = as the interface position of the steady state. 
Due to solute conservation the value of concentration on the solid side of the interface is 



Coo in the steady state, and hence Coo/ K will be the corresponding interface concentration 
in the liquid side. The steady interface temperature will be then Tm + rriLCoo/K, and the 
temperature field can be written as 

T{z) = TM + mL^ + GLZ, (18) 

where Gl is the value of the thermal gradient. 

We have finally the two boundary conditions for the concentration. The first is the solute 
conservation equation, Eq. [161 which written in dimensionless variables reads 

d,C={K-l)Cat) (19) 

The second is the local equilibrium condition at the interface (Eq. [TT]) which, after sub- 
stituting the temperature from Eq. [181 writing it at the interface {z = C(^)) and non- 
dimensionalizing takes the following form: 

{l-r)C{t) = j^{KC-l). (20) 

The parameter r is defined as: 

^ . KDGl _ Id . . 

{K-l)v,mLCoo It' ^ ^ 

where the thermal length l^ is defined as 

{K - 1) rriLCoo 
^^ = KGl • (^^^ 

It is the scale of the thermal gradient, and in our present setup is also the distance between 

the interface position at equilibrium (with Vp = 0) and the (planar) interface position in 

steady state with f p 7^ (see Fig. [1]). 

C. Time dependent solidification problem 

We will consider a linear perturbation of the flow during the initial transient of the 
directional solidification experiment. The base solution will be then a quiescent fluid, but 
with time-dependent concentration field, corresponding to the building up of the solute layer 
ahead of the planar solidification front. The final state of this transient (i.e. for t — t- 00) is 
the steady state given by 

Cst{z) = l + ^-^e-% z>0, (23) 



with the interface placed at z = 0. An exact solution for the transient without flow can 
be obtained from the numerical resolution of an integro-differential equation,— but for our 
present purposes it is more convenient to use the approximation due to Warren and Langer,- 
which is known to be remarkably accurate, ^^ and is more convenient for the stability calcu- 
lation. In this approximation we assume the ansatz that the time-dependent concentration 
has an analytical form similar to that of the steady state. We then write the concentration 
profile in the following form: 



CwL{z,t) = l + {CwL{m,t)-^)e''^, z>at), (24) 



where C\YL{C{'t),t) is the equilibrium concentration corresponding to the local temperature 
at the interface, ({t) is the interface position, and l{t) is the thickness of the solutal boundary 
layer ahead of the interface. By using this ansatz and the boundary conditions (Eqs. [T9ll20|) . 
integrating the diffusion equation for the concentration (Eq. [7]) from z = (^{t) to z = oo, 
and after some manipulation, one obtains two differential equations for the functions ({t) 
and l{tyA 



. _ 1 1 - (1 - r)C(t) 

^^'> - lit) {1 - r){K - mt) + I ^^^^ 

j(.. ^ il-r)iKat) + l{t)) _ (1 - r)l{t) 
^^ l{t) {1 + {1 - r){K - l)at)) l-(l-r)C(t) ^ ^ 



Initial conditions for these equations are ({0) = It/ Id = (1 ~ t)^^ and /(O) = 0. Eqs. 
I25ll26] are what we will call the Warren-Langer equations.- They can be seen as the lowest 
order result of an expansion in moments of the diffusion equation, and quantitatively it is 
a very good approximation^^. As we will use this solution, our calculation, being linear in 
a small perturbation of the flow, will still be zeroth order in that moment expansion. By 
numerically solving Eqs. [2511261 and introducing the results into Eq. |2ll we have an analytical 
expression for the unperturbed concentration profile during the transient that will allow for 
the instability calculation of flow in the melt. This is performed in the next section. 



III. TRANSIENT FLOW INSTABILITY 

A. Perturbative calculation 

Explicitely the perturbations for botli tlie velocity and the diffusive field are introduced 

as 

v^{x,z,t) = evi^, {z - ({t)) cos{ax)e''* (27) 

C{x, z, t) = Co{z, t)+eCi{z- C(t)) cos(ax)e'^*. (28) 

When writing the previous equations, an exponential growth for the perturbations has been 
selected, with a growth rate a. This assumption is natural, in that we will not go beyond 
linear perturbation theory. We also see that the perturbation has a wavelength a in the 
direction parallel to the interface. Furthermore, we expect the convective instability to be 
confined to the region near to the front, i.e. where solute gradient is larger, and hence 
functions Vi^z and Cq should vanish at infinity. 

It should be noted that an important approximation has been performed when specifying 
the form of the perturbations. This is the adiabatic or quasi-static approximation^ii^, and it 
has been used in many different applications^iiSI^^. This approximation consists in taking 
time as a parameter, and hence take the values of a, ({t) and l{t) as constants when solving 
the problem. The validity of the approximation is conditioned to the separation of the 
timescales associated with the (fast) growth of the perturbations and the (slow) evolution 
of the fiat interface. 

Now we introduce the previous equations into Eqs. [9] and [7] and, after linearization we 
obtain, for the linear order in e: 

((^ + «i - * (^ - "") -) (^ - "') "- - ^sj^a^C. (29) 

((1 + 0^ + ^ - "' - '^) C'l = "i^-^C-o (30) 

If we denote the differential operators on the left hand side of the previous equations by 
M. for Eq. [29] and by M for Eq. |30] we can eliminate Ci by applying M over Eq. [291 ^ind 

we obtain 

K 

UMvi^^ = a^- -Ras Vi^ACq (31) 

i — A 

10 



From f 1 2 and Eq. |29]we can recover Ci. Now, to obtain a closed solution of the previous 
equation, we substitute Co by the Warren-Langer approximation Eq. [211 Also we take ({t) 
as a known function, computable by integration of Eqs. I2M2^ Then the previous equation 
takes the form 

^fM i;i,, = -ja^Ras (1 - (1 - r)C) e'Tvi,,, (32) 

where we assume that operators and functions depend on ^ = 2; — ^(t). The operator M JH 
can be factorized in first order terms, so that the previous equation is written as 



^ ^ d ^ \ a^Ras 



A 



in which the fj take the following values (where in each case even (odd) indexes correspond 
to the plus (minus) sign choice): 

^~o,i = i ((1 + ± \/(l + C)^ + 4(a + a2)) , 

^2,3 = ^ ((1 + ± ^J{l+CY + ^Sc{a + Sca^)^ , (34) 

f4,5 = ±a. 

If we make the following changes: 

Ti = 1 - Zfj, (35) 

- ^'"'^"^(l-(l-r)C)e-f, (36) 



Sc 
and apply them successively to Eq. [331 the following equation is obtained for vi^z'- 



nB^^--0 



viA^) = 0, (37) 



which is the Generalized Hypergeometric Equation.— Solutions for that equation can be 
written with the aid of the generalized hypergeometric function, which is defined by the 
following series: 

In our case, Eq. [371 has 6 linearly independent solutions, which can be written as 

Uh{s) = s^-''\F,{- 1 + r, - r;;; s), /i = 0, 1, . . . , 5. (39) 

11 



For convenience we use a notation in which (. . . + r^ + . . .) stands for a set of arguments 
with g = 0, 1, . . . , 5. The star, e.g. (. . . + r^ — r^), indicates that the case q = h is omitted, 
then counting as 5 arguments for the function 0-^5- 

We then consider that the solution of the problem is a linear combination of the functions 
of Eqs. 123 over which we have to apply boundary conditions. First of all, boundary 
condition at infinity rules out functions Uj{s) with odd index, since these functions diverge 
when z — )■ 00, i.e. when s — t- 0. At the interface the first order of fluid velocity should verify 

vi,, = (40) 

^ = (41) 

dz 

Conditions for the concentration at the interface can be worked out to obtain 

(^^-S{t)\Mv,,, = 0, (42) 

where 

S{t) = {K-l){l + at)) . (43) 

Now, from Eqs. HQlHIl Unapplied over a linear combination of the three solution functions 
well-behaved at infinity, i.e. those with even index, we obtain a homogeneous system of three 
equations with three unknowns. By writing 

vi,z{s) = aiUo{s) + a2U2{s) + asUiis), (44) 

this system can be written as a 3 x 3 matrix C such that: 

Ci, = f/2a-i), C2,=SU2u-i), C,, = {S + lS{t))MU2U-i) (45) 

where we have used the notation S = s d/ds and M. for the operator M. once we changed the 
variables from ^ to s. In order to find a non-trivial solution, Det(G) has to be equal to zero. 
The expression of Det{C) can then be worked out through a long calculation by exploiting 
several properties of the generalized hypergeometrical functions, with the result: 

Det{C) s-3+'^o+'^2+^* = (1 _ n + l{t) S{t)) Det{A) + Det{B) (46) 

where the elements of the matrices A and B are given in the appendix \M The dispersion 
relation takes then the form 

(1 - ri + l{t) S{t)) Det{A) + Det{B) = (47) 

12 



B. Asymptotic expansion for large Rag 

The dispersion relation as derived above is of limited utility. In principle Eq. HT] can 
be solved numerically. In fact this can be done only with difficulty, and only for not too 
high values of Rayleigh number. This was already noticed by Hurle et alM in a related 
problem, who found a region in the Rag-a plane that was very hard to explore numerically. 
They called it unrepresentable singular region. Furthermore, it is difficult to extract general 
properties of these solutions, due to the fact that they depend on determinants of complicated 
hypergeometric functions. For instance, numerically one finds some modes from the multiple 
solutions of Eq. HTJ but it is difficult to generalize such result from the analytical properties 
of this equation. 

In order to obtain more useful expressions for these determinants, we perform the limit 
s — > — oo. Physically this limit corresponds to large Ras and a region not very far from the 
interface, precisely the region where the solute gradient exists and the instability is expected 
to occur. The asymptotic expansion of the generalized hypergeometrical functions in this 
limit is outlined in Appendix |Bl A cumbersome calculation provides the dispersion relation 
of the instability as the solution of the following transcendental equation: 

[1 - ri + Ns^i + lS{t)] (v^ + tan (Gs^ + 37r7o + vr (ro + rg + r^)) ) + s^ = (48) 

where A3 i has the following value: 

A3,i = ^ (2 - Sc-') (C + 1) - {^ {Sc~' (5 Sc-' - 2) + 5) (C + 1)' (49) 



- I' {Sc'' + 1) a + iy4a2+(c+l)%4a-3/V - ^ 

From Eq. |l8]we discern directly several features. On the one hand, this is a transcendental 
equation, but much simpler than Eq. HTJ and then insight on their solutions can be obtained 
by intuitively solving it for instance graphically. On the other hand the tangent is a periodic 
function. This implies directly that there will be an infinitude of solutions for this equation 
as Ras/ Sc is increased, a feature that was not immediate from the closed form of A and B 
and a property that is confirmed numerically after much difficulty for several modes. These 
solutions are bifurcations of the unstable "fiat" state, and they could have a role in the full 
non-linear dynamics of the system. 

Finally, the asymptotic expansion provides an additional result. We have proven in 
Appendix [B] that if we use only the dominant terms in the expansion, the determinants 

13 
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FIG. 2. Stability boundaries as a function of time for Sc = 81, r = 0.8, K = 0.3, and several values 
of the Rayleigh number. The dotted line is the boundary of a transiently unstable region, since its 
corresponding Rayleigh number (Ras/Sc = 9.8) is slightly below the threshold of the instability 
in the steady state {Ras/Sc = 10.28). 

cancel out exactly, and terms beyond all orders must be added to the expansion in order to 
obtain Eq. HHl This arises immediately the question of the numerical precision of the direct 
computation of the determinants of Eq. HTJ since the larger part of the numerical value of 
each matrix element will cancel out exactly. We ascertain that this is in fact the origin of 
the unrepresentable singular region of Hurle et ai— 

IV. NUMERICAL RESULTS AND DISCUSSION 

We can now solve Eq. HT] and find the stability boundary as a function of time, making 
Re{a) = (no oscillating instabilities have been found). The control parameters will be 
Ras and r. 

If we fix T and change Ras we expect that, as Ras rises, the instability will be anticipated, 
since the steady system is more unstable with higher values of Ras- Nevertheless, little can 
be said a priori about the first unstable wavelength. 
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From Fig. [2] we see a dramatic dependency of the stability balloon on Ras- For values 
slightly below the threshold, we see that a transient instability is present, and as we rise 
the Rayleigh number above the threshold this region remains detached from the steady 
instability until it eventually connects with it. 

This behavior can be easily explained if we take into account the actual form of the 
stability parameter. From the changes of variables leading to Eq. [37] we see that the 
argument of the Hypergeometric function will be indeed Ras/Sc, but supplemented with 
time-dependent factors: 

-^^^^(i-(i-.)CW) (50) 

In particular, we see that this combination of parameters depends on l{t) to the fifth 
power, hence its effect will be very important. It is known^ that values of r close to 1 
(equivalently high values of It/^d, for instance for large pulling velocities) imply that I, 
when approaching its steady state value, overshoots severely. This overshooting is therefore 
responsible for the transient instability as well as the intricate shapes of the stability balloon 
in Fig. [21 On the other hand, the value of the first unstable wavevector does not change in 
an appreciable way when we merely rise the Rayleigh number. 

We can also consider the effect of r. From the previous discussion we also expect r to be 
destabilizing, since values closer to one will mean higher overshooting. This is the case, as 
we can see in Fig. [31 We also see in this figure that, despite the first unstable wavelength 
seemed to be almost independent of the Rayleigh number, higher values of r push this 
wavevector to higher wavelengths. 

This can be understood from the values of the parameters of the hypergeometric functions 
( IMj) . Whenever a appears in the equations, it has always an / multiplying it, although not 
necessarily directly. Therefore, higher I implies smaller a. 

In principle, it might seem that a natural assumption would be to rescale all lengths with 
/ and assume that the steady states results are still valid in the transient. This was indeed 
the path followed by Jamgotchian et al.^. The Rayleigh number is rescaled by P, and a by 
/^, in accordance with the construction of the diffusion time from the diffusion length. 

We argue that this approach can only work partially. A simple calculation shows that 
when ( = the scaling is indeed perfect, but only in that case. For this value of (, 
(1 + C) = l^^, and hence the Tj in equation [351 are identical to their counterpart in the 
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FIG. 3. Stability boundaries as a function of time for Ras/Sc = 15, Sc = 81, K = 0.3, and several 
values of r. The innermost curve is the one with r = —0.1, and the values of the r parameter grow 
outwards. 



steady state. Similarly, the argument of the hypergeometric function or the product l{t) S{t) 
becomes independent of time, thus reverting to the steady-state case. 

Nevertheless, even if the transient instability is controlled solely by the value of /, the 
instability will develop in general before the time at which C = 0- I^ Fig- HI it is depicted 
the comparison of the actual time to instability, the boundary of the steady instability and 
the boundary of the transient instability as computed with the scaling approximation for 
K = 0.3 and Sc = 81. 

To begin with, we observe that for small or negative values of r the steady instability line 
is almost coincident with the boundary of the transient instability, and close to the threshold 
it takes longer for the instability to develop, as it could be expected. For larger values of r we 
see that the boundary of the transient instability goes below the critical Ras of the steady 
instability, i.e. there appears a purely transient instability, as shown in Fig. O We see that 
the scaling approximation (dashed line) is a very good one in that, despite its simplicity, 
correctly predicts the dependency of Ras with r for the instability boundary. Nevertheless, 
we see that the transient instability takes place for smaller values of Ras, as expected from 

16 




FIG. 4. Stability map for Sc = 81, K = 0.3. Color depends on the logarithm of the time to the 
first instability (blue - short, red ~ long). The solid line is the steady-state stability boundary 
of the convective instability, the dashed line corresponds to the instability when performing the 
scaling of all lengths with /. 



the previous argument that the scahng is correct for C = but only approximate when this 
is not the case. 

We can compare now the asymptotic approximation with the numerical results in the 
low a regime for some of the lowest-lying modes of the instability in the steady state (see 
Fig. |5]). We see that the asymptotic approximation correctly predicts the a~^ dependency 
for small a, and we also see that the accuracy of the approximation increases with Rag- 
Nevertheless, we see that for larger values of a the approximation is not so good. In fact, 
the asymptotic behavior for large a is not correctly reproduced. We believe that this is due 
to the fact that the arguments of the hypergeometric functions, which in turn depend on 
the Ti, are the same in the large a limit. This corresponds to a degenerate problem, and our 
method cannot be applied in such circumstances. 

Finally, the consistency of the approximations performed should be tested against the 
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FIG. 5. Successive instabilities of the homogeneous state (modes 2, 3, 4, 5 and 6) for Sc = 81, 
K = 0.3. Sohd: Numerical solution. Dashed: Asymptotic approximation. 



numerical results. In particular, it should be checked that the adiabatic or quasi-static ap- 
proximation is indeed a good one. Obviously, near the threshold of the instability this is 
not the case, since at that point the evolution of the perturbation is extremely slow (a ~ 0). 
Therefore, the location of the threshold will not be precisely determined, but the approxi- 
mation will still be useful to locate regions where perturbations will be either amplified or 
damped. To test the validity of the approximation and, in turn, discuss the observability 
of the transient instability, we have computed the growth rate a of the perturbation for the 
most unstable mode in the transient instability found for Rag/ Sc = 9.8 and r = 0.8 (see 
Fig. 12]). The value found is arD = 0.14 for t/ro = 10. Hence, even in such an extreme 
situation, the separation of scales is valid {t ^ l/o")- In addition, at that point the derivative 
of a with respect to time is close to zero, being as it is close to a maximum. Therefore, the 
time evolution of the fiat interface is slow enough for the adiabatic approximation to give 
consistent results and for the transient instability to develop. 
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V. CONCLUSIONS 

We have studied the convective stabihty of the melt during the initial transient of a 
directional solidification experiment in a vertical configuration. The setup is such that 
thermal gradient is stabilizing but the solute layer built during the transient is destabilizing. 

By using a quasi-stationary approximation we have obtained the time-dependent disper- 
sion relation for perturbations (both in fluid velocity and in concentration) acting over the 
quiescent solution in which the solute field follows the Warren-Langer approximation for the 
transient of a planar solidification front. Such dispersion relation has a very complicated 
analytical structure and appears as very hard to evaluate numerically. We have further per- 
formed an asymptotic expansion for large Rayleigh numbers which has permitted on the one 
hand to extract some properties of the solution and on the other hand to explain the origin 
of the difficulties that affected the numerical calculations. In particular the asymptotic pro- 
cedure reveals that the dominant terms of the dispersion relation cancel out exactly. The 
result comes from sub-dominant terms, and hence is much harder to evaluate numerically 
from the complete solution. The analytical form of the asymptotic solution also shows the 
existence of an infinite number of bifurcations of the unstable state. 

The numerical evaluation of the dispersion relation has permitted to work out the stability 
boundary on the parameter space as a function of time. As it could be expected large values 
of r (i.e. small thermal gradients) tend to destabilize the fiow, and to displace the instability 
to larger wavelengths. In this regime, increasing Rayleigh number to above the threshold 
a transient instability appears, i.e. the quiescent solution is stable for the steady (large 
times) state, but it has not been so during the transient. This instability could deform the 
solidification front and drive the system to a very different final state, which would not have 
been revealed in a purely steady state analysis. By increasing Ras an instability at large 
times appears, and by increasing it further both instability regimes merge. This could be 
relevant for microstructure predictions in solidification processes. 

In this study we have not considered the possible deformation of the interface. In fact 
the solidification front does eventually undergo a morphological instability, which in the 
directional solidification setup occurs for velocities above a threshold value.- In this insta- 
bility the transient considered here turns out to be of capital importance. ^^ An interesting 
continuation of this work would be the complete analysis of the coupling between the mor- 
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phological and the convective dynamics (already analyzed for the steady state2^>^) during 
the transient. Such analysis would permit to discern the regimes in which both dynamics 
effectively couple to each other in the initial transient, presumably giving rise to different 
results from what given by each of them separately. 
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Appendix A: Explicit expressions for matrices A and B 

The elements of the matrices A and B can be written as: 

All =Bii= oF^i; 1 + ri - ro, 1 + r2 - ro, 1 + rs - tq, 1 + r4 - tq, 1 + rs - tq; s) 
Ai2 =Bi2= oi^5(; 1 + ro - ra, 1 + ri - r2, 1 + rg - rs, 1 + r4 - rg, 1 + rg - rg; s) 
Ai3 =Bi3= oF^i; 1 + ro - r4, 1 + ri - r4, 1 + ra - r4, 1 + r^ - r^,l + r^ - r^; s) 

A21 = B21 = (ri - ro)oF5(; ri - ro, 1 + rg - ro, 1 + rg - ro, 1 + r4 - ro, 1 + rg - ro; s) 
A22 = B22 = (ri - r2)oF5(; 1 + ro - r2, ri - r2, 1 + rg - r2, 1 + r4 - r2, 1 + rg - r2; s) 
A23 = B23 = (ri - r4)oF5(; 1 + ro - r4, ri - r4, 1 + r2 - r4, 1 + rs - r4, 1 + r5 - r4; s) 

5 
A31 = n('"i ~ ro)oF^^{; 1 + ri - ro, r2 - ro, rg - ro, r4 - ro, rg - ro; s) 

i=2 
A23 = 

(l + ro-r4)(l + ri - r4) 

5 
B31 = n(^i ~ ^o)o^5(; n - ro, r2 - ro, rg - ro, r4 - ro, rg - ro; s) 

_ s oF5(; 2 + ro - r2, 1 + ri - r2, 1 + rg - r2, 1 + r4 - r2, 1 + rs - r2; s) 

^33 — 7:r- ^ 

(1 + ro - r2) 

So-F5(; 2 + ro - r4, 1 + ri - r4, 1+ r2 - r4, 1 + rs - r4, 1 + rs - r4; s) 



SoF5(;2 + ro- 


- r2, 2 + ri - r2, 1 + rs - r2, 1 + r4 - 


- r2, 1 + rs - 


-r2;s) 


SoF5(;2 + ro- 


(1 +ro-r2)(l + ri - r2) 
- r4, 2 + ri - r4, 1 + r2 - r4, 1 + rs - 


- r4, 1 + rs - 


- ^4; s) 



B 



33 



(1 + ro - r4) 

Note that although the elements of A and B are defined as functions of s they are to be 
evaluated ai z = ({t) 
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Appendix B: Asymptotic expansion of the dispersion relation 

Conventional asymptotics for the pFq function (p 7^ 0) for large arguments makes use 
of the Mellin-Barnes representation. By shifting the integration contour, the sum of the 
residues on the poles of T{aj + s) gives an asymptotic expansion for large argument^, 
ignoring exponential terms beyond all orders.— When j9 = 0, there are not that kind of 
poles, and hence the expansion is necessarily exponential. ^° By working out this expansion, 
particularized for g = 5, we can reach the following expression, in which the expansion can 
be written as the sum of three series: 



with 



nr(P. 



.1=1 



4V3^ 



(5o + 5i + ^2) , (Bl) 



5o(p,; z) = z^e'^'^'^' ^ ^^^^^^ ^^ " Y + ^''^J ""'' ^^2) 



00 o , ^^ 



Slip,; z) = -z^T.12 Nrcos 37r7 - — + 62^/^ + 2npj ^'i, (B3) 

r=Oj=l V / / 

S^ip,; z) = z^e-'^^'^' £ J2 iV.cos (bTT^ - ^ + Sz'^' + 27c{p, + pk)) z'l . (B4) 
In these series, the coefficients Nk fulfill the following recursion formula 



s— 1 


A;>0 


iVo = l, iV, = 0, 


A; <0; 



kNk = Y.T^-s{s-k)Q~''^Nk-s, k>0 (B5) 

(B6) 



where 



Tit) = l[(i + ^3) , (B8) 

a;, = 6(p,-l + 7). (B9) 
We follow the notation in which p^ = 1, and the parameter 7 is given by 

T-^-igp.. (BIO) 

We next use this expansion to find an approximation in the limit s — )■ — 00 of the deter- 
minants of the matrices A and B. To do so, some shorthand notation is needed. We will 
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consider a general matrix whose elements are generalized hypergeometric functions. If the 
element {i,j) of one of these matrices is o-Fg(;Pij,i, ■■■,Pij,5', s) then the quantities defined in 
Eqs. IB6|B7|B8|B9I can be computed for each (i, j), by taking these specific values of Pg, and 
hence can be re-defined as matrices. 

Now, when introducing these series into the matrix elements of A and B given in Ap- 
pendix IX], it can be seen that most of the factors accompanying the generalized hypergeo- 
metric functions can be absorbed, after some manipulations with the series, into the factor 
11^=1 r(pj), which becomes the same for all the elements. There are still a —z factor on some 
elements of the third row. This can be eliminated by redefining their specific jij values with 
the substitution 732 — )■ 732 + 1, 733 — ;■ 733 + 1, where the prime denotes the new value of 
gamma. With these substitutions a simple calculation shows that, for matrix A: 

72, = 7u + g, (Bll) 

73, = 71, + I- (B12) 
For matrix B: 

72, = 71, + I, (B13) 

73, = 71, + g- (B14) 

The result is that, discarding the global factor 11^=1 r(Pi), every element (i, j) of both ma- 
trices A and B can be written as a sum of three contributions Sij^, Sij^i, Sij^2, each one 
written as in Eqs. IB2|B3IIB4t but with the understanding of taking the element {i,j) of each 
quantity depending on the Pg. 

The resulting series depend on the N,. coefficients, which theirselves depend on the actual 
parameters of the hypergeometric function through a recurrence relation which involves a 
polynomial whose roots are Uj = /3(pj — 1 + 7). We are now going to use the actual values 
of the Pg to prove that in fact for each of the functions the roots are the same up to a 
permutation. 

To do that we can compute the matrix elements of Uk for each determinant and each row. 
First of all, particularizing for the first row of both A and B matrices, it can be seen that the 
set of the {uJk}k=o,i...5 are in fact just permutations for different j. Therefore, since the roots 
of the polynomial can be exchanged, the value of N^^ does not depend on j. For the second 
row of both A and B we notice that the coefficients of the hypergeometric functions on that 

22 



row are in fact the same if we make the substitution ri -^ r[ + 1. We can then proceed as 
in the previous row and prove that N2j,r does not depend on j either, and its value is the 
same. Finally, for the third row, we have to distinguish between A and B. For matrix A, the 
substitutions Tq — t- Tq — 1, ri — )■ r'^^ — 1 will bring us to the first row, and on matrix B, this 
will be accomplished by the substitution rg — > Tq — 1. Therefore A^3j,r is independent of j, 
and the value is also the same as the other rows. The final conclusion is that the coefficients 
of the asymptotic expansion A^^ do not depend on the particular eigenfunction. 

We are now ready for expanding the determinants of the matrices A and B. We introduce 

the following notation: 

Sij,i{r) 

S2j,Us) (B15) 

where Sijj{r) stands for the whole row, not just the j^^ element. Written in that way we 
can express the whole series, by using the properties of the determinants, as 



'j-rst 
I h 



Imn 



|A| -EEEEEEt;: 

I m n r s 



rst 
Imn 



(B16) 



1. Dominant terms 



With the previous definitions, the dominant terms can be denoted by T^oo* • '^^^ corre- 
sponding determinant has the form: 

f + 3^1/6^ cos (7172 -f + 3^1/6) cos (7r73 -f + Sz^^^) 

(vTTi - ^^ + 3^1/6) cos (7r72 - ^^ + 3z'/^) cos (7r73 - ^^ + 3z'/^) (BIT) 



K 



cos [ TTJi 
COS 



COS 



(^71 - 6"^ + 3^-*^/^ j COS f 7r72 



7r(f-4) 



+ 



3zi/6) COS (^73 - ^^ + 3zi/6j 



where 



K = N^,rN2,sN^,tZ 



7i+72+73-(r-+s+t-5)/6^9V32l/6 



(B18) 



We have used the definitions of the 7jj for A (this choice is immaterial, as it will be clearly 
shown) and the substitutions A^j ,, = Ajj^r, Ij = lij- Expanding the cosine functions of the 
second and third rows, and using the properties of the determinants, it is very simple to 
show that 

7^00 = (B19) 

i.e., all the dominant terms in the asymptotic expansion of the determinant cancel exactly. 
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2. Beyond all orders 

Since all the dominant terms cancel, sub-exponential terms are needed. The first sub- 
exponential terms are of the form: 

V^o (B20) 

These terms will be computed in their general form, but before some more notation is 
required. We make the following definitions: 

70 = -^-^Er. (B21) 

tto = -^ (B22) 

o 

«i = ^ (B23) 

4 - 1 

a2 = (matrix A) (B24) 

5 — t 

as = (matrix B) (B25) 

6 

K = Ni^rN2,sN3^tZ^'+^'^''''-^'+'-^'-^'^/^e^^''^' (matrix A) (B26) 

K = Ni^rN2,sN3,tz'"^^'+^''+^-^''+'+'^/^e^^''^' (matrix B) (B27) 

With these definitions, the general determinant can be written as: 

5 

Tloo = -Keijk J2 ^^^ (S^^ + 3vr(7o + ao) + 7rr2(i_i) + 27rn_e(2(i-i)-o) (B28) 

1=1 

cos (Sze + 7r(7o + tti) + 7rr2(j_i)j cos (Sz'^ + 7r(7o + 0:2) + 7rr2(A:_i)) , (B29) 



in which summation over i,j,k from 1 to 3 is implied, eijk is the Levi-Civita symbol, and 
6{n) is a discrete step function with 6{n) = for n < and 6{n) = 1 for n > 0. 

After a lengthly manipulation of the trigonometric functions, and by using symmetry 
properties under exchanging of indexes, one can simplify this last expression to get 

7^00 = -4-ft'sin (7r(ai - 02)) sin (tt (r2 - r^)) sin (vr (r4 - tq)) sin (tt (ro - r2)) 

cos [Qz'^ + 37r7o + 37rQ;o + tt (^0 + ^2 + ^4) ) • (B30) 



Notice that in fact we have only computed a separate term. We have assumed a certain 
position for the sub-exponential term in the determinant, but to compute properly the 
contribution we should sum for all possibilities: 



'100 ^ '010 ^ '001 ' 



7-rst I '-j-rst I '-j-r 
inn T /Qin ~r 'nr 
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This is tantamount to a cyclic sum over the ctj. Further manipulation gives 



T'rst _i_ '-i-rst _i_ '-i-rst 
inn -r /qio "T '001 



'100 



— IG/Tsin {n(ai — 0:2)) sin (7r(a2 — Q^o)) sin (7r(a;o — en)) 

sin (tt (r2 — Vi)) sin (tt (r4 — tq)) sin (tt (ro — r2)) (B31) 

cos (^62:6 + Svr^g + vr (tto + ai + ^2) + tt (ro + r2 + r4) 



We can now evaluate the terms corresponding to the first orders in the large- 2; expansion. 
For each matrix A and B we consider first the zeroth order terms: 



-7-000 I -7-000 I -7-000 

nnn "r /mn "r /qoi • 



'100 



'010 



Note that those terms are not of the same order in z for the matrices A and B, due to the 
different order of the K for each matrix. Matrix B gives in fact the lower order in z. 

We start then with the matrix B. For this matrix the value of the Oj are a 
which in Eq. IB31I gives: 



^66 



Tir + ToT + To^ = 2v^sin (tt (r2 - r,)) sin (vr {r, - ro)) sin (vr (ro - r2)) 

cos (Gz^ + 37170 + TT (ro + r2 + r^)) ;27i +72+73+1^6^3.1/6 ^g^^) 



Order zero of matrix A contributes to the next order in z. For this matrix the value of 



the Oj are given by a 



012 



. If we put that in Eq. IB31I we obtain: 



7^00° + 'Cio" + T^oT = -4y3sin (tt (r2 - n)) sin (tt {n - ro)) sin (vr (ro - r2)) 

^l + ^(r,^r.4- r.)\ ^7i+72+73+5/6,6V3.V6 
6 



cos ( Qz^ + 37170 + ^ + TT (ro + r2 + r4) ) 27i+72+73+&/egt>V3.-" (333) 



We only consider now the next order of matrix B, which is of the same order in z than 
the zeroth order of matrix A. There are nine of those order 1 terms: 



-7-100 
'100 



-7-100 
'010 



-7-100 
'001 



-7-010 I -7-010 I -7-010 I -7-001 
/inn -r /mn "t" /ooi "t" 'lOO 



'100 



'010 



-7-001 
'010 



-7-001 
'001 



We can evaluate the previous expression by summing terms in groups of 3. For 7^oo° + 
Tqiq + Tqqi the vector of the a^ is a = — | | | • Since 0^2 — Oq = 1 one of the sines in Eq. 
IB31l will be zero and this term will not contribute. The next group is 7^oo° + '^010° + '^i°- 



Here, the vector of the a,- is a 



00 



. For this case ao ^ c^i = and again one of the 



sinus in Eq. IB31l will be zero and this term will neither contribute. For the group of terms 
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7?m + T^m + T£?? the vector of the a^ is a 



'100 ^ '010 ~ '001 



12 
"63 



. But this is the same vector that was 



previously found for A at order 0. The final result will then be, for the matrix B: 
r,Z + 'T^io' + T^oi' = -4V3iV3,isin (vr {r^ - n)) sin (vr {n - ro)) sin (vr {r^ - r,)) 

cos Uz^ + 37r7o + ^ + ^(^0 + ^2+ n)) 2^i+^2+73+5/6g6v^.i/6 

with A^3^i as in Eq. HH 

Finally, we have all the elements at our disposal to build an approximation to the dis- 
persion relation in the form given by Eq. |46l If we substitute our approximation in this 
formula we obtain: 

- 4^3 [1 - ri + Ns^i + ISit)] cos (6zh3n-fo + y + vr (ro + rg + n)) (B34) 

+ z'^2v3cos(Qz'^ + 3iT'yo + Tc {ro + r2 + r^)] = 
After rearranging the terms we can write it in a more compact form: 

[1 - ri + A^3,i + lS{t)] (VS + tan {Gz'^ + 87170 + tt (ro + r2 + n))) + z'^ = (B35) 

Which, clearly, is a transcendental equation for z. The direct analysis of equation |48] shows 
that there exists an infinite number of solution branches as z increases, due to the periodicity 
of the tangent. 
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